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Abstract 


This  report  describes  NRL’s  hydrodynamic  (isopycnal)  nonlinear,  primitive 
equation,  ocean  circulation  model  in  a  spherical  layer  suitable  for  a  large  basin 
or  the  global  domain.  This  document  should  be  read  in  conjimction  with 
NOARL  Report' 35,  December  1991,  The  Navy  Layered  Ocean  Model  User's 
Guide  [25]. 

Spherical  geometry  requires  careful  treatment  of  the  equations  of  fluid  dy¬ 
namics  to  ensure  spurious  terms  are  not  introduced  during  the  layer  averaging 
process,  particularly  in  the  momentum  difiusion  expressions. 

Prom  these  layer  averaged  equations,  finite  difiierence  approximations  are 
derived  on  a  grid  of  uniform  intervals  in  longitude  and  latitude  (not  necessarily 
the  same).  These  equations  may  be  solved  efficiently  on  modern  architecture 
scientific  computers,  allowing  accurate  description  of  ocean  flows. 
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IV 


EXECUTIVE  SUMMARY 


This  report  describes  the  formulation  of  NRL’s  hydrodynamic  (isopycnal)  nonlinear,  primitive 
equation,  ocean  circulation  model  in  a  spherical  layer  suitable  for  a  large  basin  or  the  global  domain. 
This  document  should  be  read  in  conjunction  with  NOARL  Report  35,  December  1991,  The  Navy 
Layered  Ocean  Model  User's  Guide  (hereafter  NOARL35)  [25]. 

Spherical  geometry  requires  careful  treatment  of  the  equations  of  fluid  dynamics  to  ensure 
spurious  terms  are  not  introduced  during  the  layer  averaging  process,  particularly  in  the  momentum 
difliision  expressions. 

From  these  layer  average  equations,  flnite  difference  approximations  axe  derived  on  a  grid  of 
uniform  intervals  in  longitude  and  latitude  (not  necessarily  the  same).  The  resulting  flnite  differ¬ 
ence  equations  may  be  solved  efficiently  on  modern  architecture  scientiflc  computers  to  describe 
accmately  ocean  currents,  temperature  and  salinity  distributions  over  large  time  and  space  scales. 


E-1 


Formulation  of  the  NRL  Layered  Ocean  Model 
in  Spherical  Coordinates 


> 


1.  PRELIMINARIES 

1.1  Motivation 

Ocean  Modeling  has  evolved  in  20  years  from  simple  limited  resolution,  limited  area,  ideal 
domain,  Cartesian  geometry  computations  [7],  [21],  [8]  to  realistic  basin  scale  eddy  resolving  ocean 
models  [6],  [10],  [9].  This  progress  has  reflected  both  the  growth  in  computing  capacity  and  in 
the  sophistication  and  experience  of  numerical  oceanographers.  However,  the  actual  models  cannot 
form  a  continuous  range  as  the  transition  from  Cartesian  to  spherical  geometry  results  in  substantial 
changes  in  the  appearance  of  the  vector  calculus  operators  and  of  the  equations  of  fluid  flow,  in 
spite  of  the  coordinate  invariance  of  the  underlying  physics  and  mathematics.  This  change  in 
appearance  reflects  the  difficulties  of  expressing  the  laws  of  physics  in  a  cmrvilinear  coordinate 
system,  particularly  the  conservation  of  vector  quantities  such  as  momentum. 

Another  difficulty  is  that  while  the  choice  and  orientation  of  {x,y,z}  to  label  the  orthogonal 
Cartesian  coordinates  appears  to  be  common  across  most  physical  disciplines,  no  such  consensus 
exists  for  the  choice  and  orientation  of  spherical  or  spheroidal  coordinates.  While  most  mathemat¬ 
ical  and  theoretical  studies  favour  co-latitude  as  the  declination  coordinate,  many  oceanographers 
and  meteorologist  prefer  latitude.  However,  introductory  mathematics  and  physics  text  books 
tend  to  opt  for  the  former  choice.  This  can  lead  to  confusion  and  error.  The  labels  assigned  to  the 
thre4  spherical  surface  coordinates  also  vary  widely,  even  within  the  meteorology  and  oceanography 
literature. 

For  simplicity,  the  criteria  of  minimizing  differences  between  the  Cartesian  and  spherical  models 
will  guide  the  choices  of  notation  and  orientation  used  in  this  presentation. 

A  note  of  caution.  Care  will  need  to  be  taken  to  distinguish  two-dimensional  mathematical 
operations  defined  on  the  surface  of  a  sphere  and  three-dimensional  physical  and  mathematical 
processes  defined  within  a  thin  spherical  layer.  The  purpose  of  the  Shallow  Water  set  of  approxi¬ 
mations  is  to  reduce  the  equations  that  apply  to  three-dimensional  physical  flow  throughout  a  thin 
fluid  layer  to  equations  describing  a  two-dimensional  model  flow  on  a  spherical  surface  in  a  way 
that  preserves  as  much  of  the  physics  of  the  3-D  flow  as  possible. 

1.2  Coordinates  and  unit  vectors. 

To  conform  to  the  conventional  right-handed  ix,y,z)  Cartesian  coordinate  system  locally,  we 
chose  longitude,  latitude  and  radius  as  our  surface  spherical  coordinate  directions  in  that  order. 
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We  label  the  longitude  and  latitude  angles  (f>  and  6  respectively.  We  assume  that  <j)  increases  in 
an  eastward  direction,  and  that  6  increases  in  a  northward  direction  from  a  value  of  zero  at  the 
equator,  taking  negative  values  in  the  Southern  hemisphere.  The  two  spherical  surface  coordinate 
singularities  (the  North  and  South  Poles)  are  assumed  to  coincide  with  the  rotation  axis  of  the  earth, 
unless  stated  otherwise.  The  origin  of  the  longitude  (^)  coordinate  can  be  chosen  for  convenience, 
with  the  Greenwich  meridian  assumed  to  be  the  default.  The  surface  of  the  earth  is  chosen  to  lie 
on  a  sphere  at  a  constant  radius  a.  The  radial  coordinate  is  chosen  to  increase  outwards.  We  will 
label  this  coordinate  r.  As  the  angles  (f>  and  6  are  dimensionless,  physical  distances  will  be  given 
implicitly  in  the  form  acf)  or  a  0  ,  where  a  is  a  suitable  mean  reference  radius  for  the  earth. 

We  label  the  unit  vectors  that  axe  locally  tangent  to  the  surface  and  parallel  to  lines  of  constant 
latitude  and  constant  longitude  as  and  Sq  respectively.  We  assume  that  points  eastward 
and  bq  points  northward.  There  is  a  local  unit  normal  vector  pointing  out  of  the  sphere  at  each 
point  (0, 0,  a)  on  the  spherical  surface. 

The  distinction  between  a  sphere  and  true  geopotential  surfaces  will  be  ignored.  Gill  [4]  pp  91  - 
92  contains  a  full  discussion  of  the  error  incurred  in  approximating  the  shape  of  the  earth  by  a 
sphere.  There  he  estimates  that  a  maximum  error  of  0.17%  occurs  in  the  geometrical  factors  in  the 
various  equations  if  the  earth’s  figure  is  approximated  as  a  sphere  instead  of  an  oblate  spheroid. 
Thus,  although  in  reality  horizontal  fiuid  fiows  are  locally  parallel  to  geopotential  surfaces,  very 
little  error  is  made  approximating  these  complicated  surfaces  by  concentric  spheres. 

1.3  History 

The  derivation  of  the  spherical  Shallow  Water  equations  from  the  full  3-D  Navier  Stokes  equa¬ 
tions  in  a  spherical  coordinate  system  has  been  carried  out  by  several  authors  in  the  past  [11], 
[17].  While  early  agreement  was  reached  as  to  the  proper  form  for  the  momentum  advection  and 
conservation  of  mass  expressions,  opinions  continue  to  differ  as  to  the  correct  form  for  the  layer 
averaged  viscous  dissipation  in  the  presence  of  layer  thickness  and  horizontal  density  variations. 
We  followed  Wajsowicz  [24]  and  Shchepetkin  O’Brien  [18]  in  the  derivation  of  these  terms  and 
can  demonstrate  that  our  forms  generate  no  spurious  forces  for  analytic  test  cases  such  as  rigid 
rotation  of  a  varying  thickness  shell. 

Given  the  complexity  of  the  appearance  of  these  equations  in  spherical  coordinates,  many 
mathematically  equivalent  expressions  are  possible.  These  variations  in  appearance  can  lead  to 
variations  in  the  finite  difference  approximations  used  when  replacing  these  partial  differential 
equations  by  difference  equations.  These  can  give  rise  to  widely  differing  answers  to  the  same 
underlying  mathematical  problem.  We  note  that  while  the  equations  given  in  Section  2  are  mathe¬ 
matically  equivalent  to  those  appearing  in  Shchepetkin  and  O’Brien  (1996)  [18],  the  most  recently 
published  contribution  to  this  discussion  that  we  are  aware  of;  they  are  not  identical.  They  differ 
in  the  form  taken  for  the  viscous  dissipation,  although  they  adhere  to  the  requirement  that  this 
term  be  the  divergence  of  a  symmetric  tensor.  The  form  we  use  is  taken  from  Wajsowicz  (1994) 
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2.  SPHERICAL  SHELL  FLUID  DYNAMICS 

The  reduction  of  the  three  dimensional  partial  differential  equations  describing  the  conservation 
of  mass,  momentum,  energy  and  matter  (tracers)  at  each  point  in  an  ocean  layer  to  a  set  of  two- 
dimensional  partial  differential  equations  for  flow  on  a  spherical  layer  is  done  by  ntegrating  these 
equations  in  the  radial  direction  over  the  thickness  of  the  spherical  shell,  h  from  r  =  a  tor  =  a  +  h. 
This  averaging  process  requires  application  of  boundary  conditions  at  the  top  and  bottom  of  each 
layer  and  the  conversion  of  surface  stresses  into  body  forces.  In  our  equations  these  are  subsmned 
into  the  “somces  —  sinks  ”  terms. 


In  Cartesian  coordinates,  the  integration  is  straightforward  and  yields  the  equations  set  out 
in  NOARL35  [25]  or  Hurlburt  and  Thompson  [8].  On  a  spherical  shell  the  averaging  must  be 
done  with  care  to  avoid  generating  spurious  sources  or  sinks  of  energy,  mass  or  momentmn  from 
curvature  terms.  After  a  lively  discussion  in  the  literature  (sf.  [1],  [2],  [5],  [11],  [12],  [13],  [15]  [17] 
[20]  [23]  &  [26])  a  consensus  has  emerged  as  to  the  correct  forms  for  these  equations  and  we  have 
chosen  a  representation  that  closely  resembles  the  one  used  in  our  earlier  Cartesian  work  [25]. 

2.1  2-D  Spherical  Shell  Velocity  Field 


Classically,  the  derivation  of  the  Shallow  Water  equations  assumes  each  layer  is  homogeneous 
in  the  vertical  direction  and  that  vertical  velocities  within  the  layer  can  be  ignored  except  to 
balance  horizontal  flow  divergence  and  conserve  mass.  While  this  vertical  uniformity  assmnption 
is  appropriate  in  Cartesian  coordinates,  it  introduces  spurious  shear  forces  in  a  spherical  layer.  It 
is  necessary  to  assume  that  the  horizontal  velocity  has  a  radial  dependence.  This  ensures  that 
particles  at  the  top  and  bottom  of  the  layer  that  lie  on  a  line  through  the  center  of  the  sphere 
have  the  same  angular  velocity  about  the  center  of  the  sphere  and  hence  remain  in  step  for  any 
horizontal  velocity  field. 


The  appropriate  form  for  an  incompressible  velocity  field  (u)  in  a  thin  spherical  layer  is: 


where  A  &  B  are  chosen  to  match  the  conditions: 

w{4),0,r  =  a,t)  =  wi,{(l),d,t), 
w{(f>,e,r  =  a  +  h,t)  =  Wb{<f>,e,t)  -I- 

at 


(2a) 

(2b) 


This  linear  radial  dependence  of  the  horizontal  components  of  the  velocity  field  is  necessary  to  keep 
the  top  of  a  column  of  fluid  positioned  over  its  bottom  for  all  horizontal  flows.  The  entire  layer 
moves  coherently  as  it  flows  over  the  surface  of  the  sphere  and  there  is  no  shear  between  fluid  at 
different  radii.  The  vertical  velocity  is  chosen  to  match  the  vertical  motion  of  the  lower  boundary 
and  the  local  divergence  of  the  horizontal  motion  field.  This  maintains  the  incompressibility  of  the 
overall  flow. 


For  reference: 

w{<f>,9,r,  t) 


(r^  -  Q^)  (a  +  hf%  -f  [(r^  -f  ha^)  (a  +  hf  -  r^a^]  Wb 
r2[(a  +  h)3  -  a3] 


(3) 
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In  the  limit  <g  1: 


Ar^+B 

J.2 


Wb{i>,6,t)  + 


di’ 


similar  to  the  Cartesian  result  [25]. 


The  leading  error  term  is  : 


/rj-a\  /a  +  h  —  r\  r dh  /h 
\  h  )  \  a  )  _  dt  ^  Va 


Wb 


This  vanishes  at  r  =  a  and  r  =  a  +  h  and  is  bounded  in  that  interval  by; 


1 

4 


(4) 


(5) 

(6) 


The  divergence  of  this  velocity  field  (l)is: 


OCOS0 


du'  d(v'  cos  9) 


3{a  +  h)^  dh  ^  3  [{a  +  h)^  -  g^] 


(o  +  h)^  —  dt  {a  +  h)^  — 


Wb. 


(7) 


As  this  relationship  does  not  depend  upon  the  vertical  coordinate,  it  is  valid  everywhere  in  the 
fiuid  layer. 


In  the  limit  0,  we  recover  (for  an  incompressible  fluid)  a  two-dimensional  continuity 

equation  : 


or 


—  +  1  ( d{hu')  d{hv'  cos0)\  ^ 

dt  a  cos  6  \  d(j)  89  }  ^ 


(8) 


dh 

—  -h  V2-F  «  0,  (9) 

(using  the  subscript  2  to  denote  spherical  surface  two  dimensional  vector  calculus  operators)  where 
V  is  the  horizontal  velocity  transport  vector  (see  next  page ). 


The  leading  error  term  is 


(10) 


when  approximating  (7)  by  (9).  In  the  absence  of  vertical  motion  of  the  lower  boundary,  (9)  is 
accurate  to  O  . 
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2.2  A  Single  Layer  with  Uniform  Density 

Here  are  the  equations  used  within  a  single  layer.  If  this  is  one  of  a  stack  of  layers,  interactions 
between  layers  are  mediated  by  the  “sources  -  sinks  ”  terms  in  each  equation.  The  effective  pressure 
terms  ^  will  have  to  be  modified  as  in  NOARL35  [25]  to  reflect  the  possibilties  of  barotropic 
and  baroclinic  modes  of  motion.  Hereafter  we  will  drop  the  prime  to  distinguish  2-D  velocities  from 
3-D  velocities  and  assume  that  (u,  v)  label  horizontal  layer  averaged  velocity  components. 

2.2.1  Continuity 


dt  a  cos  6 


m  djV  cos  9) 
dcj)  86 


=  sources  -  sinks. 


(11) 


2.2.2  Velocity  Transport 

Define  an  angular  deformation  tensor  e  to  have  components: 


8  1 

(  u  \ 

f.  8  1 

'  ^  ^ 

(  =  -  egg) 

^tfxp  - 

8(1)  ' 

^COS0/ 

,COS0/  ’ 

(12a) 

5 

/  «  \ 

.  S  1 

^  u  \ 

( =  eg^) 

e<t>e  = 

50 

\cos0y 

<COS0/  * 

(12b) 

where  the  velocities  {u,  v) 

are  defined 

from  the 

transports  {U,  V)  by: 

hu  =  U, 

(13a) 

hv  =  V. 

(13b) 

Then  the  two  components  of  the  Velocity  Transport  equation  may  be  written: 
component: 


dU  1 
-t- 


dt  a  cos  6 


diUu)  ,  diVucosO)  _  .  ^ 

— 50 - ^  - 50 - —  aflV  sin20 


hg  dh 


+  hF^  +  -2^2 

^  cos^  9 


d{h  cos  6  d{h  cos^  9  e^g) 

d(t>  5F^ 


a  cos  9  50 
+  sources  -  sinks,  (14a) 


eg  component: 

dV  1 

+ 


dt  a  cos  6 

Ah 


d{Uv)  ,  5(Vt;cos0)  ^  ' 

— 50 ^  50 - U u  sm 9  +  aQ,U  sin 20 


hg  dh 
a  89 


-\-hFg  + 


COs2  0 


8{h  cos  0  e^)  d{h  cos^  0 
50  50 


-h  sources  -  sinks. 


(14b) 
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2.2.3  Scalar  transport 


d{hT)  1  \d{UT)  diVTcosd^ 
dt  "^acos^L  d<f>  d9  . 


K  r_^  / 

0?  cos^  6  ,d(j>\  d<f>) 


Q 

+  COS«^ 


+  sources  -  sinks. 


(15) 


Here  h  is  the  vertical  layer  thickness;  f2  is  the  rotation  rate  of  the  earth;  g  is  the  acceleration  due 
to  gravity;  and  Fg  axe  the  components  of  the  applied  forcing  fields  in  local  (^,  6)  components; 
Aff  is  the  eddy  viscosity;  T  is  a  scalar  tracer;  and  k  is  the  coefficient  of  diffusivity  for  the  scalar  T. 


2.3  Discussion 


These  forms  are  mathematically  equivalent  to  those  used  in  Shchepetkin  and  O’Brien  (1996) 
[18].  The  form  of  the  diffusion  operator  is  taken  from  Wajsowicz  (1994)  [24],  assuming  that  the 
vertical  diffusivity  and  second  viscosity  coefficients  vanish.  This  form  of  the  friction  term  should 
be  more  accurate  for  layers  with  varying  thickness  than  the  form  used  in  some  earlier  Cartesian 
studies;  in  Hmlburt  and  Thompson  (1980)  [8],  although  it  will  require  more  work  to 

calculate. 


Equations  (12a-b)  imply  a  stress-strain  relationship  of  the  form 


,  O'  =  Ah  h  cos  9  e, 

and  a  stress  force  contribution  to  the  transport  equations  of  the  form: 

1 


a  cos  9 


V2  a, 


(16) 


where  V2‘  is  defined  in  Appendix  B  (B5).  cr  and  e  are  angular  tensors.  These  require  a  factor 

of  3  to  transform  them  to  physical  tensors  [22].  The  components  of  I  may  be  related  to  the 
— » 

rate-of-strain  tensor  £  components  (as  defined  in  [3])  for  the  velocity  field  (1): 

^<p<t>  —  0'{£(pif>  Sgg) ,  epg  =  1a£^.  (17) 

Equations  (11)  -  (15)  are  valid  only  in  the  limits  <  1  and  «,  u  »  ty.  Typical  terms  ignored 
when  deriving  these  simplified  equations  are  of  the  order  of  f-)^. 
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2.4  Enstrophy  Conserving  Formulation 

An  alternative  formulation  (as  recommended  by  [17],  [1])  uses  the  radial  vorticity  component 
C  defined  as  : 


.  _  1  r  ^  5(ucos0)'l 

a  cos  9  \d(j)  do  /  ’ 

and  exploits  the  vector  identity 

(u  •  V)«  =  (V  X  n)  X  u  +  V  . 


This  substitution  yields: 

dU  1 
“Trr  + 


dt  a  cos  9 


.  dU  ^  dV  ,  d(Vcos9) 


-  (C  +  2fisin^)  V, 


(18) 


(19) 


(20a) 


^  1 
dt  ^  a 


dV  dU  V  tsaOVv 
2v  —  +  u  —  + 


do  do  cos  0  a 

as  an  alternative  left  hand  side  for  equations  (14a)  and  (14b). 


+  (C  +  2fl  sin  0)  U, 


(20b) 


If  u  is  used  instead  of  F  as  the  prognostic  variable,  then  the  left  hand  sides  of  these  equations 
assume  the  simpler  forms  [Ij: 


h  d  f  +  v‘^ 


du 


dt  a  cos  0d4>  \  2 


-  (2ftsin0  +  0^ 


,  dv  h  d  ( +  v‘^\  .  _ 


(21a) 


(21b) 


i 
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2.5  Relocating  the  Poles 

Converging  coordinate  lines  may  cause  difficulties  with  many  numerical  schemes,  particularly 
in  limiting  the  size  of  the  permissible  time  steps.  A  ploy  to  lessen  this  problem  is  to  move  the 

coordinate  singularities  so  that  they  both  lie  over  land  [19].  If  this  is  done,  the  Coriolis  term  in  the 

momentum  transport  equations  needs  to  reflect  this  loss  of  coincidence  of  the  rotation  axis  and  the 
North  and  South  Poles  of  the  coordinate  system. 

In  vector^forn^  the  Coriolis  pseudo-force  in  the  velocity  transport  equations  (14a,  14b)  is  ex¬ 
pressed  as  2  il  X  U .  Using 

ft  =  (cos^q  sin0q  i  +  sin^Q  sin0q  j  -h  sin^q  k) 

and  re-expressing  ft  in  the  local  spherical  basis  fusing  equations  (Ala-c)  from  Appendix  A)  we  find; 

ft  =  ft  cos^q  sin(^q  -  4>) 

+  ft(sin0q  cos<^  —  cos^q  sin0  cos(^q  —  (f))  eg 
4-  ft  (sin^q  sin0  -t-  COS0q  COS^  cos  (^q  -  (f)  )  Cr, 

=  (ft(^,  fig,  fir) . 

For  [/  =  (U,  V,  0)  the  Coriolis  pseudo-force  vector  becomes  2  (-Uftr,  Ufir,  Yfi^f,  -  Ufig) 
or,  in  spherical  surface  coordinates,  {—V  f,  U  f)  where: 

/  =  2ft  (sin0q  sind  -|-  cos^q  cos0  cos(0q  —  (p)) .  (22) 

Thus,  on  a  shifted  coordinate  spherical  surface,  the  Coriolis  term  becomes  dependent  on  both  the 
local  longitude  and  latitude,  but  retains  its  2ft  (— U  /,  U  f)  form. 
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3.  FINITE  DIFFERENCE  APPROXIMATIONS  TO  THE  EQUATIONS 
3.1  The  grid  on  a  sphere 


Finite  difference  approximations  to  the  equations  set  out  in  Section  2  are  defined  on  a  uniform 
(A^,  A6)  grid  covering  most  of  the  earth,  but  commonly  excluding  extreme  latitudes.  The  ’C’ 
grid  of  [13]  defines  the  spatial  staggered  grids  where  the  h,  U  and  V  variables  take  discrete  values. 
Figure  1  sets  out  the  indexing  of  these  variables  on  the  three  interlocked  staggered  meshes. 
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Figure  1 

The  relative  layout  of  the  h,  U  and  V  grids. 

The  dashes  mark  the  boundary  of  the  fluid. 

The  normal  components  of  velocity  vanish  across  this  line. 

A  fourth  grid,  vorticity  (Cij)j  can  be  defined  at  the  points  (02i+i ,  92j+i ),  lying  horizontally  between 
the  Vij  points  and  vertically  between  the  Uij  points.  ^  ^ 
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S.1.1  Infinitesmal  vs.  finite  area  elements  and  line  lengths 

An  infinitesmal  difierential  surface  area  element  on  the  sphere  in  our  coordinates  is 
a^  cos  9  d9d^.  The  true  area  of  a  finite  mesh  cell  centered  at  {<p,e)  of  width  a  cosOAcf),  and 
height  aA$  is  2a^  A<^  sin  cos  6  .  To  leading  order,  the  proportional  error  between  the  true 

area  and  the  finite  area  element  c?  cos  6A(f)  A9  is: 

true  area  -  finite  element  A6^ 
true  area  ~  24  " 


More  elaborate  finite  diflference  schemes  than  those  set  out  in  the  subsequent  subsections  of 
this  paper  that  intend  to  achieve  higher  accuracy  will  have  to  incorporate  this  second  order  finite 
area  difference. 


Fortunately,  finite  line  elements  parallel  to  the  (<^,  9)  coordinate  axes  remain  analogs  of  the 
infinitesmal  line  elements;  (a  cos  9  dcj),  a  d9)  becoming  exactly  (o  cos  9  A<^,  a  A9)  . 

3,1.2  Point  vs.  Area  or  line  averaged  values 


The  variable  is  a  layer  thickness  value  centered  at  the  point  +  i  A(f> , 

=  Oq  +  jA9  averaged  over  a  mesh  cell  of  width  acos9jA(f>  and  height  aA9  .  The  variable 
is  an  average  of  a  transport  over  a  line  from  the  point  {(f>i,9j)  to  the  point  i(f>i,9j+i)  . 
Similarly,  is  defined  as  the  average  of  the  transport  across  the  line  from  {<pi,9j)  to  the 

point  dj)  .  In  our  spherical  coordinate  system  these  averages  differ  from  spot  values  centered 

on  {(i>i,9j)  by  terms  of  second  order  and  higher  in  A9,A<p  times  <p  and  9  derivatives  of  the 
field. 


If  we  consider  in  detail  the  relationship  between  a  spot  value  of  a  field  /i(^o,^o)  and  the 
average  value  of  this  field  over  the  finite  area  element  of  size  2a'^  A(f>  sin  cos  9  centered  at 
the  point  (^O)  &o)  {called  h  )  we  find: 


h  «  hi<f,o,9o)  + 


^  (A0)2 

'd'^h 

r.  .  dh 

■ 

+  — — 
.0  24 

d9^ 

-  2  tan^o 

So 

(t>0  ^ 

+  0{A(t>\A<l>'^A9'^,A9^)  + 


(24) 


A  similar  second  order  dependence  on  the  mesh  spacing  exists  between  the  average  value  of 
U  between  the  points  {(f>if9j)  and  {<pi,  0j+i)  and  its  spot  value  U{<l)i,  9j_^x)  •  Likewise  for  V  , 

Again  these  differences  between  spot  and  average  values  must  be  refiected  in  more  accurate 
difference  formulae. 
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3.2  Second  Order  Finite  Difference  Equations 

For  convenience  in  notation,  let  us  follow  other  authors  in  ocean  modelling  [11]  and  define 
centered  averaging  and  difference  operators: 


mA\ 

-w(z-  ^))  /(mA), 

(25a) 

mA\ 

(25b) 

The  centered  time  t  applies  to  all  fields  unless  otherwise  noted.  Second  order  accurate  finite 
difference  approximations  to  the  equations  set  out  in  Section  2.2  take  the  following  forms: 

3.2.1  Conservation  of  Mass 

Equation  (11): 


52th  = 


-1 

a  cos  6 


S^U 


+  Sg{V  COS  9) 


4-  sources  —  sinks. 


(26) 


3.2.2  Velocities 


The  velocities  {u,  v)  are  defined  from  the  transports  {U,  V)  by: 


4 

3.2.3  Deformation  Tensor  Components 


(27) 


Equations  (12a)  and  (12b): 


(28a) 

(28b) 


This  defines  the  diagonal  elements  of  the  deformation  tensor  and  egg  on  the  hij  grid  and 
the  off-diagonal  elements  and  on  the  Qj  grid.  This  positions  these  values  exactly  where 
they  are  needed  by  the  difference  operators  to  calculate  the  viscous  stress  contribution  to  the  U 
and  V  advection  equations. 
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3.8.4  Velocity  transport 

Equations  (14a)  and  (14b)  are  approximated  by: 

^2tU  =  +  S0{y^xf  cosO^  —  sinO  {if  h)'^ u  —  a Q  sin  26  {if  h)^ 

Q  ~h^  d)  U  [  1  ^  ^ 

^cos^g  +  Se(lif  008“^ 6 e^')  ,  (29a) 


S2tV  =  6tj,(u^lf’'^  +  6g(y^lf  cos^j  +  {sinOUu^ h)  +  afl (sin  2du’l‘hf 

—  ^~^^9h+  fi  Fq  +  (;Qg2  Q  cos  0  —  5e^/icos^0e<^^jj  .  (29b) 

3.8.5  Scalar  tracers 

Scalar  fields  such  as  temperature,  salinity  and  tracer  concentrations  (T,  S',  Ci)  are  defined  on 
the  hij  grid.  Equation  (15)  may  be  approximated  in  a  manner  consistent  with  the  previous  finite 
difference  equations  by: 

S2t{hT)  =  ^64UT’^)  +  Se{VT^  cosO) 

<  +  (.Qg2  0  +  COS  9  5$  (h  COS  e^SeT)]'  (30) 

3.8.6  Discussion 

The  horizontal  diffusion  terms  in  the  three  finite  difference  equations  given  above  are  lagged  in 
time  one  time  step  for  computational  stability  [16].  This  introduces  first  order  timestepping  errors. 
In  practice,  small  values  of  v  and  k  are  used,  in  fact  the  smallest  that  allow  the  code  to  be  stable. 
Hence  a  small  constant  proceeds  this  first  order  timestep  error  and  its  overall  effect  is  hoped  to  be 
small. 

In  the  momentum  equations  for  time  stepping  U  &  V,  the  advection  terms  have  a  conserva¬ 
tive  flux  contribution  and  separate  curvature  terms  (—V u  tand,  U  u  tan0).  These  curvature  terms 
redistribute  the  momentum,  but  if  the  equation  {u  x  (14a)  -|-  v  x  (14b))  is  formed  to  pro¬ 
duce  a  prognostic  equation  for  a  kinetic  energy  (u  17  -f  u  F),  it  is  seen  that  these  curvature  terms 
make  no  contribution  to  the  kinetic  energy  evolution,  just  as  the  Coriolis  terms  also  cancel  in  this 
combination. 
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Holland  and  Lin  [5]  devised  finite  difference  forms  for  the  Coriolis  terms  in  each 
equation  {-{v^fh)  ,{u'^fh)  )  that  preserve  this  conservation  of  kinetic  energy  when 
—  {U^fhfv^  is  summed  over  a  domain  where  the  normal  velocity  component  van¬ 
ishes  at  the  boundary.  This  may  be  confirmed  numerically  for  arbitrary  (/,  h,  17,  V)  fields  satisfying 
the  boundary  conditions.  When  the  axis  of  rotation  goes  through  the  coordinate  singularities 
(North  &  South  Poles),  this  finite  difference  approximation  reduces  to  the  form  set  out  above. 
The  more  general  forms  given  in  this  paragraph  should  be  used  when  the  poles  are  moved.  Other 
suggested  finite  difference  approximations  for  the  Coriolis  terms  [18]  do  not  conserve  kinetic  energy. 

This  energy  conservative  finite  difference  approximation  to  the  Coriolis  terms  is  chosen  as  a 
guide  to  approximating  the  curvature  term  in  each  equation. 

The  scalar  tracer  finite  difference  equation  preserves  locally  the  property  of  no  change  in  the 
scalar  field  when  it  is  constant,  no  matter  what  the  velocity  field.  Any  local  divergence  in  the  17,  V 
field  is  refiected  only  in  a  change  of  layer  thickness. 

3.3  Diagnostic  Quantities  and  Energetics 

3.3.1  Vorticity 


The  normal  component  of  vorticity  is  defined  on  a  spherical  surface  grid  [4]  as: 

.  _  1  t  dv  d{u  cos  6)  I 

a  cos  6  \  d(j)  dO  /  ’ 


(31) 


The  finite  difference  analogue  of  this  equation  is: 


Cm  =  ~  Seiucosd)}. 


It  may  be  used  as  a  diagnostic  variable. 

3,3,2  Kinetic  Energy 

If  the  kinetic  energy  field  is  defined  as  UijUij  +  VijVij,  the  two  products  are  defined  on 
separate  grids.  One  or  both  need  to  be  extrapolated  onto  a  common  grid  to  define  a  kinetic  energy 
field,  although  the  separate  grids  may  be  used  for  the  purposes  of  integrating  kinetic  energy  over 
the  entire  domain.  If  the  kinetic  energy  is  to  be  combined  or  compared  with  the  potential  energy, 
then  the  definition  on  the  hij  grid  would  be  appropriate,  viz,: 

iKE)ij  = 


(33) 
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4.  TRUNCATION  ERROR  AND  ACCURACY 

As  the  operators  (25a, 25b)  used  in  constructing  the  finite  difference  approximations  to  the 
equations  of  fiuid  motion  are  centered,  we  expect  them  to  represent  these  differential  equations 
with  error  terms  proportional  to  even  powers  of  the  mesh  spacings  (A^,  A0)  and  their  products. 
A  full  truncation  error  analysis  of  the  finite  difference  equations  in  section  3.2  will  reveal  this 
expected  behavior. 

It  is  useful  to  attempt  to  quantify  these  error  expressions  by  assuming  a  surface  fiow  of  a 
simple  type  and  seeing  how  these  finite  difference  equations  reproduce  the  properties  of  the  analytic 
equations  for  this  fiow.  The  explicit  form  taken  by  the  truncation  errors  in  this  case  may  be  a  useful 
guide  as  to  the  choice  of  mesh  spacing  and  to  the  degree  the  distortion  of  the  mesh  near  the  poles 
degrades  the  accuracy  of  these  equations.  One  such  simple  fiow  is  that  of  a  rigid  shell  in  unform 
rotation  about  some  pole  centered  at  the  point  {<I)r,  Or)  . 

4.1  Analytic 

The  surface  velocity  field  (u,  v)  of  a  thin  spherical  shell  rotating  from  west  to  east  about  an 
axis  passing  through  the  point  ((pR,  Or)  with  angular  velocity  u)  is: 

Uij  =  oa;[cos(^  -  (^r)  cos  Or  sin^  -sin^jR  cos^,  sin(<^H  -  cos^ij].  (34) 


If  UR  is  used  to  prescribe  17  and  F  in  the  continuity  equation  (11),  it  can  be  shown  that  the 
expression  in  the  square  brackets  vanishes  identically.  If  ffij  is  inserted  into  the  expressions  for  the 
viscous  stress  tensor  components  then  it  be  shown  that  they  all  vanish  everywhere,  as  expected  for 
such  a  fiow. 

If  these  forms  are  used  for  17,  V,  u  and  v  in  the  bilinear  advection  terms  for  the  velocity 
transport  time  evolution,  the  terms  that  arise  can  be  balanced  exactly  by  the  pseudo-Coriolis  force 
for  an  axis  of  rotation  moved  to  (  4>Ri  Or)  and  a  rotation  rate  of  —  a; .  This  is  as  expected  when 
the  apparent  forces  arising  from  a  rotation  of  the  inertial  axes  is  cancelled  by  an  apparent  fiow  field 
in  the  moving  frame  exactly  opposite  to  the  rotation  field. 

These  analytic  results  will  be  useful  for  investigating  the  accuracy  of  finite  difference  represen¬ 
tations  of  the  governing  equations  and  alternative  formulations. 
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4.2  Numerical 

If  the  above  velocity  field  is  projected  onto  a  grid  of  size  (A^,  A^)  and  these  discrete  values 
are  used  in  the  finite  difference  approximations  for  the  governing  equations  we  can  look  at  the 
differences  between  the  discrete  results  and  the  analytic  results  in  our  several  cases. 

4-2.1  Continuity 


The  right  hand  side  of  equation  (26)  evaluates  to: 


cos  Or  sin^  tan^ 


(35) 


Thus  if  A^  —  A0 ,  the  terms  inside  the  square  bracket  cancel  and  our  finite  difiierence  scheme  is 
exact  for  this  model  fiow.  The  difference  between  A^  and  A^  determines  the  actual  trimcation 
error  for  this  scheme  and  grid.  To  the  lowest  order,  this  error  takes  the  form: 


52th 


cos  Or  sin  (f)  tan  6 


(A^)" 

24 


^  +  o((A^)4-(A0)4)  +  ... 


4-2,2  Advection  terms 


If  the  general  rigid  rotation  velocity  field  (34)  is  substituted  into  the  advection  part  of  the 
transport  equations,  very  cpmplicated  expressions  result.  The  general  level  of  accuracy  of  the  finite 
difference  approximations  here  can  be  illustrated  by  considering  two  limits  of  the  rotation  velocity 
field;  (i)  axis  of  rotation  at  the  poles  and  (ii)  axis  of  rotation  at  the  equator. 

I^et  the  velocity  field  for  case  (i)  be  u  =  (-  cosO,  0) .  This  velocity  field  in  the  co-rotating 
frame  exactly  recovers  the  inertial  frame  and  there  are  no  real  forces  in  the  advection  terms.  By 
inspection  we  can  see  that  the  finite  difference  approximations  to  the  advection  terms  in  the  U 
transport  equation  (29a)  evaluate  to  zero  for  any  choice  of  A<p  and  AO.  In  the  V  tranport 
approximation  (29b)  the  curvature  term  is  cancelled  exactly  for  any  mesh  by  the  Coriolis  force 
finite  difierence  approximation. 


^  When  we  move  the  axis  of  rotation  to  the  equator  at  =  0 ,  the  velocity  field  takes  the 
form:  u  =  ( —  cos sin0,  sin^) .  The  analytic  contribution  to  the  U  transport  advection  is 
sin^  cos  <p  cos^  0 .  The  finite  difference  approximation  evaluates  to: 


sin^  cos(f)  cos^  9 


,  3  +  14tan^^  . 

1 - 24 - ^ 


7  cos  20 
24  cos^  0 


AO^  + 


(36) 


Both  of  the  coefficients  of  the  second  order  terms  are  well  behaved  between  70°  N  and  70°  S ,  being 
bounded  by  3.3  for  the  first  coefficient  and  -1.25  for  the  second  in  this  range.  The  truncation  error 
properties  of  the  finite  difference  approximation  to  the  V  transport  equation  (29b)  are  similar. 
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4-2,3  Deformation  tensor 

iProm  the  components  (25a-b)  for  the  deformation  tensor,  we  recover  truncation  error  expres¬ 
sions  of  the  form 

^(t>(i>{^R)  «  cos^ij  sin0  tan^ 


e^5(uii)  «  cosOr  cos(f>sece  +  (2  +  6  tan2(0)^  +  ...  .  (38) 

Although  the  factors  involving  6  tan^^  can  get  large  at  high  latitudes,  these  approximations  are 
seen  to  retain  very  good  second  order  accuracy  over  much  of  the  globe. 

When  these  expressions  are  used  in  calculating  the  stress  contribution  to  the  velocity  treins- 
port  equations  (26a-b),  the  resultant  truncation  errors  retain  their  ((A(^)2,  (A0)2)  character  and 
magnitude. 
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Appendix  A 

THE  CONNECTION  BETWEEN  CARTESIAN  AND  SPHERICAL  BASIS 

VECTORS 


It  is  useful  to  be  able  to  re-express  vectors  defined  in  an  absolute  (FJ,  k)  Cartesian  basis  in 
terms  of  the  local  spherical  basis  ( e^,  eg,  er)  at  the  point  {(f),  6,  r)  and  vice  versa.  The  nine 
components  of  the  matrix  that  transforms  a  vector  representation  from  one  basis  to  the  other  can 
be  expressed  in  either  the  local  spherical  coordinates  {(j),  6,  r)  or  in  the  local  cartesian  coordinates 
(x,  y,  z).  We  will  assume  for  convenience  that  the  origins  of  the  two  coordinate  systems  coincide 
and  that  the  two  spherical  coordinate  singularities  (the  North  and  South  Poles)  lie  on  the  ^-axis. 

A.l  Spherical  coordinate  representations 

To  move  firom  a  Cartesian  basis  to  a  spherical  basis  we  use  the  relations: 

60  =  —  sin  <f>  i  -|-  cos  (j)  j 

60  =  -cos0sin0  r  -  sin^sin^  j  +  cos 9  k, 

Or  =  cos^cosd  r  +  sin^cos^  J  +  sinO  k. 

To  move  fi:om  a  representation  in  the  spherical  basis  to  a  representation  in  the  Cartesian  basis 
we  use  the  relations: 

4 


• 

1 

=  —  sin^  60  —  cos^sin^ 

ee 

+ 

cos  ({>  cos  6 

©r 

(A2a) 

J 

=  cos^  60  —  sin  ^  sin  0 

60 

+ 

sin  <!>  cos  6 

©r 

(A2b) 

k 

=  COS0 

+ 

sin^ 

gr. 

(A2c) 

If  these  transformation  rules  are  expressed  as  Vj  =  T  vq,  Vg  is  the  3  components  of  v  expressed 
in. the  local  spherical  basis  set  (  e^,  e^,  Cr),  vc  is  that  same  vector  expressed  in  the  Cartesian 

basis  set  ^i,  J,  k  j  and  T  is  a  three  by  three  matrix.  We  note  that  ^T^  =  ^T^  .  The  matrix 

T  is  orthogonal  as  its  transpose  is  its  inverse! 


(Ala) 

(Alb) 

(Ale) 
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A.2  The  Spherical-Cartesian  transformations  in  Cartesian  Coordinates 
— ♦ 

To  express  T  purely  in  terms  of  local  Cartesian  variables,  we  define  r  =  , 

p  =  (x^  +  Then 


X 

P 

r  p 


\ 


X 

T 


u 

r 


(A3) 


In  this  form  it  is  easy  to  demonstrate  that  the  basis  (e^,  eg,  e^)  defined  by  (Al  a-c)  has  the 
(required  !)  properties: 


II 

1, 

II 

0, 

*  ©r 

=  0, 

= 

^6'  Gg  = 

1, 

er 

=  0, 

(A4) 

X 

II 

Sr  X  eg  = 

*  Gr 

=  1. 

Appendix  B 

SPHERICAL  SHELL  VECTOR  CALCULUS 


B.l  Definition  of  Surface  Differential  Operators  on  a  Sphere 
B.l.a  Basic  Operators 

Let  $  be  any  scalar  quantity,  and  let  it  be  a  function  of  (^,  $)  only, 

(viz.  —  =  0). 

Let  F  be  a  vector  quantity  such  as  velocity  v,  or  mass  transport  V.  F  is  assumed  to  lie  in  the 
local  tangent  plane  and  to  describe  a  field  everywhere  nearly  parallel  to  the  spherical  layer  it  is 
embedded  in.  Quantities  advected  by  v  (or  V)  are  assumed  to  maintain  vertical  coherence  across 
a  fiuid  layer  (c.  f.  [14]  p  63).  That  is,  two  particles  lying  on  the  same  radial  vector  at  the  top  and 
bottom  of  a  fluid  layer  respectively  remain  on  a  single  radial  vector  as  they  are  advected  by  the 
velocity  field  of  that  layer.  This  implies  that  u  is  obliged  to  assume  the  form: 

u(<f),  e,r)=(^^  e),  ^  v'{(f),  6),  w(<f>,  e,r)^,  (Bla) 

where  w((f>,6,r)  takes  the  form  necessary  to  maintain  zero  divergence  for  the  velocity  and  to 

match  the  motions  of  the  upper  and  lower  boundaries  of  the  layer,  in  what  follows,  let  us  consider 
a  general  vector  F  is  expressed  in  coordinate  form  as: 

F  =  ^F^{<f>,e),Fr(<t>,e,r)^ .  (Bib) 

In  many  ocean  model  contexts,  the  radial  component  of  the  velocity  vector  is  ignored.  For  these 
purposes  it  is  useful  to  consider  just  the  horizontal'  components  of  velocity.  Let  us  define  U2  to 
be: 

v'{4>,e),  0) ,  (B2) 

and  similarly  F2  to  be  any  locally  horizontal  vector  field. 

The  standard  invariant  operators  of  vector  calculus  are  now  set  out  in  explicit  coordinate  form 
using  our  spherical  surface  variables  (<^,  0)  and  assuming  a  functional  form  for  vectors  defined  as 
in  (Bib)  or  (B2)  and  for  scalars  dependent  upon  (^,0)  only.^ 


^c.f.  Batchelor  (1970),  pp  598-601  or  Gill  (1982),  pp  92-93 
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B.l.b  Gradient  (V2$) 


V2$ 


1  ^ 

rcosO  r  89 


(B3) 


This  operator  retains  an  r  dependence  throughout  a  spherical  shell  and  will  have 
averaged  with  care. 

B.l.c  Divergence  (V  •  F) 

V-F  =  ‘  gWcoaO)  1  d(r^Fr) 

a  COS  9  d(f>  a  cos  9  89  r"^  8t  ' 


to  be  vertically 


(B4) 


The  first  two  terms  are  independent  of  r.  This  expression  may  be  used  to  define  an  r  component 
to  a  vector  field  to  make  it  divergence  firee.  We  will  define  V2'  to  be: 


a  cos  9  8<f>  a  cos  9  89 


(B5) 


B.l.d  Curl  (V  X  F) 


4 


W  X  F 


^  _ L_ 

Vr  50  a  )  y  o  r  cos 9  8<l>  j 

1  f  ^  _  d{P^cos9)  1  ^ 

a  cos  9  \  8((>  89  j 


ee 


(B6) 


When  F2  is  the  velocity  field,  the  is  the  normal  component  of  the  spherical  surface  fiow  vorticity, 
commonly  labeled  C-  Note  that  there  axe  both  and  eg  components  of  the  full  vorticity  vector 
for  all  non-zero  horizontal  velocity  fields,  even  in  the  absence  of  any  vertical  velocity. 


B.l.e  Scalar  Laplacian  (V^^) 


vi$  = 


1  8H  1  5/5$  \ 

cos^  0  cos  6  89  \  36  )  * 


(B7) 


^Pedlosky  (1987),  p.  55,  eqn.  (2.9.20)  is  incorrect.  Consider  Stokes’s  theorem  : 

( V  Xu)*  ds  =  u  •  for  the  infinitesmal  surface  area  element  A(j>  cos  9  A9  . 
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B.l.f  Vector  Laplacian  (V^^) 

1 


= 


+ 


ra  cos^  B 

1 

racos^B 

2 

ra  cos  B 


■d% 

„  d 

dcfP' 

+  cose^ 

(““^  m  ) 

'd^F^ 

Q 

+  cosB—  1 

fcosS^) 

d(jP 

dB  ' 

1  dB  J 

9F^  d{FgCosB)'\  ^ 
d<t>  dB  ®'“ 


(B8) 


We  recognize  the  er  component  as  —  ^  Va  •  F2. 

B.2  Higher  Vector  Derivatives  of  the  velocity  field 


B.2.a 


Momentum  Advection 
{U2  •  V)  U2  = 


U 


tand 


a 


u'v' 


r 

+  - 


dv!  v'  du' 
a  Locos0  d(f)  ^  a  dB 

^  t  i 

o  Locos 0  5^  a  dB  ^  a 


ee 


«  o  +  vv 


a 


B.2ib 


curl  curl 


V  X  ( V  X  U2)  =  — 


1 


+ 


ra  cos^  B 
1 

ra  cos^  B 
2 

ra  cos  B 


0  +  2coa^«„' -  |(cos«^) 


(B9) 


—  COS0 


d^v' 


d<f>dB 


6^ 


d{v'  cos  B) 
d<j)  dB 


er. 


(BIO) 


Note:  V  X  (V  X  U2)  -^^2^2  because  ^2-^2  #  0. 

B.2.C  grad  div 

,  ^  \  ^  (  1  I  du'  ^  5(w'cos0)  1 ^ 
ra[dB  \cos0  I  5^  J )\ 

+  0  •  er 


(Bll) 
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B.2.d  div  grad 


(V-V)«2  = 


1  [ d‘^u'  ,  (  ^du> 

ra  cos^  6  dcjr  dO  \  dO 

+  0  •  er 


(B12) 


B.3  The  rate-of-strain  tensor 


Assuming  w  is  not  zero: 


^<p(f}  ^tp9  ^(jyr 

Se<i>  See  Ser 

Sr(f>  Sj-e  Srr 


(B13) 


Err  — 


(B14) 


c  _  c  _  dw 
Or  -  Ere  -  2r  dd' 


(B15) 


Assuming  w  is  zero: 


p  _  p  _  •*•  'JUJ 

04^  -  Cr4  - 


1  dw 
2rcos0  d<f)' 


P  1  dv' 


(B16) 


(B17) 


Edd,  = 


1  du'  v'  tan  6 
acosO  dcp  a  ’ 


(B18) 


Sipe  =  £94  = 


cos  6  /  u'  \  1 

d$  Vcos0/  2a  cc 


1  dv' 
cos  6  d4>' 


(B19) 


